#include "HY1C_out.h"
/* ========================================================================
 * MScaltarget - finds L2 sensor overlap over a L1B scene, generates 
 *		L2 target and aerosol binaries 
 * 
 * Synopsis:
 *
 *   MScaltarget par=parfile
 *
 * Description:
 * 
 * Modification history:
 *
 *     Programmer     Organization      Date       Description of change
 *   --------------   ------------    --------     ---------------------
 *   Ewa Kwiatkowska   SAIC        24 July 2003    Original development
 *   Ewa Kwiatkowska   SAIC        24 July 2005    Matching L3s with a L1B scene
 *   Ewa Kwiatkowska   SAIC      6 November 2006   New SeaWiFS to MODIS nLw conversion
 *   B. Franz          NASA      26 September 2011 Added climatology test
 *                                                 masked more unused code
 *
 * ======================================================================== */

#include <stdio.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <stdlib.h>
#include <fcntl.h>
#include <string.h>
#include <unistd.h>
#include <libgen.h>
#include <math.h>
#include <float.h>
#include <time.h>

#include "l12_proto.h"
#include "seabin.h"
#include "readL2scan.h"
#include "l12_parms.h"
#include "l2_flags.h"
#include "hdf.h"
#include "target_struc.h"
#include "l1_struc.h"
#include "aer_struc.h"
#include "filehandle.h"
#include "input_struc.h"
#include <setupflags.h>

// int32 read_attrs(int32 sdfid, meta_l3bType *meta_l3b);

#define  RADIUS_MULTPL   5
#define  OUTSCANLINE     100  /* writing out 100 scan-lines if cannot allocate memory to the entire scene */
#define  MAXNVDATA 	 50


int32 close_input_l3b(int32 fid, int32 sdfid, int32 vgid, int32 vdata_id[MAXNVDATA]);
int32 open_input_l3b(char *hdf_file, int32 *fid, int32 *sdfid, int32 *vgid, int32 vdata_id[MAXNVDATA], meta_l3bType *meta_l3b, int32 *nrows);
int   windex(float wave, float twave[], int ntwave);
typedef int32   l3nvdata[MAXNVDATA];


int32_t  diffday(int32_t date1, int32_t date2);
void  meanday(int32_t year1, int32_t day1, int32_t msec1, int32_t year2, int32_t day2, int32_t msec2, int32_t *myear, int32_t *mday, int32_t *mmsec);
uint64_t diffmsec(int32_t year1, int32_t day1, int32_t msec1, int32_t year2, int32_t day2, int32_t msec2);
int   nLw_SeaWiFS2MODIS(int32_t sensorID, float *nLwbands, float *Lwbands);
void  tau_interp(int32_t nl1b, float *l1wave, int32_t nl2b, float *l2wave, float *l2tau);
int   compar(const void *x, const void *y);
int32 binlocate(int32 *binlist_buf, int32 n_records, int32 bin);
int   read9km_mask(char *file, char *mask);
int   is_masked(int32_t bin, char *mask, int32_t nrows);
void  bandpass_correct(int32_t sensorID1, int32_t sensorID2, float *wave1, float *wave2, int32_t nwave1, int32_t nwave2, float *nLw, float *Lw);

int empty_tgtaer( int32_t, int32_t, int32_t, char *, char * );



/* -------------------------------------------------------------------- *
 *                              main                                    *
 * -------------------------------------------------------------------- */
int main(int argc, char *argv[])
{

    int32_t   iscan    = 0;          /* input scan number                  */
    int32_t   oscan    = 0;          /* output scan number                 */
    int32_t   spix     = 0;          /* start pixel for subscene process   */
    int32_t   epix     = -1;         /* end pixel for subscene process     */
    int32_t   ipix;
    int32_t   sscan    = 0;          /* start scan for subscene process    */
    int32_t   escan    = -1;         /* end scan for subscene process      */
    int32_t   opix     = 0;          /* output pixels number               */
    int32_t   npix     = 0;          /* Number of output pixels per scan   */
    int32_t   nscan    = 0;          /* Number of output scans             */

    int32_t   climatology = 0;

    int32_t   ldate[3], mdate[2][3], fpixs;
    uint64_t  udate[3], *timediff, *xtimediff;
    char   il3file[FILENAME_MAX];

    l1str       l1rec;              /* generic level-1b scan structure      */
    filehandle  l1file;             /* input file handle                    */
    int32_t        nl2files, tl2files; /* number of l2 calibrator files        */
    int32_t        l2hdr[5]; 
    int32 	*numbin;
    int32 	*basebin;
    int32 	bin_row, bin_col, bin, ibin, bin_num;
    int32 	nrows, flags;
    int32 	n_records, interlace, interlace0, vdata_size, vsize;
    
    int32_t 	sd_id;
    int32 	*fid;
    int32 	*sdfid;
    int32 	*vgid;
    l3nvdata 	*vdata_id;
    int32 	binindex_buf;
    meta_l3bType   *meta_l3b;
    float32 	*sum_buf[2], ***sum_prod, **sum_chl, f32, v32, wgt, wgtsum;
    uint8   	*binlist_buf=NULL;
    int32   	**bins, *nbins;
    
    int16   	nobs;
    int32_t 	l2sensorID, l1sensorID, *l1wave, *l2wave, ib;
   
    char        l2sensor[400];
    int32_t        prodinx[400], chlinx, aotinx, pcvalid, badconv;
    int32_t        nprods, nprod, nl2bands, nl1bands, nl2vbands, nl1vbands, *inx;
    int         *l2finx, status;
    
    float       minLon, maxLon, minLat, maxLat, minLat1, maxLat1, lat1, lon1;
    float 	s, *nodata, chl;
    int32_t 	*d2s;
    int32_t tglength, arlength, ll, ff, scan100, wscan;
    unsigned char  *tout, *aout;
    tgstr       tgtstr;
    aestr       aerstr;
    
    float   lat, lon, *l2wavef, *l1wavef;
    int     bndflg, bndflg1;
    int     lonTest, close2thepole, west2east;
    int     latTest, close2thepole1, west2east1;

    uint32_t *flagusemask;
    uint32  required_mask;
    char    buf[5000], *mask, maskname[1000];
    char    data_prodname[128];
    char    fields[60], vdata_name[H4_MAX_NC_NAME];
    char    *outtarget, *outaer;
    instr   input;
    FILE    *outtgt, *outar;
    
    int32    i, j, n, k, l, no, format, interp_swf_mod;
    
    time_t  initialTime;
    time_t  finalTime;

    HY1C_out("\n");


    if (argc == 1) {
      l2gen_usage("vcaltarget");
      return 1;
    }

    // see if help on command line
    for(i=0; i<argc; i++) {
        if((strcmp(argv[i], "-h") == 0) || (strcmp(argv[i], "-help") == 0)) {
            l2gen_usage("vcaltarget");
            return 1;
        }
    }

    filehandle_init(&l1file);
    msl12_input_init(&input);

    /* Parse input parameters */
    if (msl12_input(argc, argv, "vcaltarget", &input, &l1file) != 0) {
        HY1C_out("-E- %s: Error parsing input parameters, exiting...\n",argv[0]);
        exit(FATAL_ERROR);
    }
    
    if (input.ifile[0][0] == '\0') {
        HY1C_out("-E- %s: No L1 file provided, exiting...\n",argv[0]);
        exit(FATAL_ERROR);
    }


    /* Determine number of input L2 files. */
    for (i=0, nl2files=0; i<MAX_OFILES; i++, nl2files++)
      if (input.il2file[i][0] == '\0')
        break;

    tl2files = nl2files;

    if (input.tgtfile[0] == '\0') {
        HY1C_out("-W- %s: No name given to L2 target binary file\n",argv[0]);
	strcpy(input.tgtfile, input.il2file[0]);
	strcat(input.tgtfile, "_Lw.bin\x0");
	HY1C_out("        Assuming name %s", input.tgtfile);
    }
    outtarget = (char *)input.tgtfile;
    
    if (input.aerfile[0] == '\0') {
	strcpy(input.aerfile, input.il2file[0]);
	strcat(input.aerfile, "_aer.bin\x0");
    }
    outaer = (char *)input.aerfile;

    if (access(input.ifile[0], F_OK) || access(input.ifile[0], R_OK)) {
        HY1C_out("-E- %s: Input file '%s' does not exist or cannot open.\n",
               argv[0], input.ifile[0]);
        exit(FATAL_ERROR);
    }
    
    for (i=0; i<nl2files; i++) {
    	if (access(input.il2file[i], F_OK) || access(input.il2file[i], R_OK)) {
        	HY1C_out("-E- %s: Input file '%s' does not exist or cannot open.\n",
               		argv[0], input.il2file[i]);
        	exit(FATAL_ERROR);
	}
    }
    

    /*									*/
    /*									*/
    /* Open input file and get sensor and scan information from handle  */
    /*                                                                  */
    format = l1file.format;
    if (format == FMT_L1XCAL) {
        HY1C_out("-E- %s: Incorrect input file format.\n",argv[0]);
        exit(FATAL_ERROR);
    }
    if (openl1(&l1file) != 0) {
        HY1C_out("-E- %s: Error opening %s for reading.\n",argv[0],l1file.name);
        exit(FATAL_ERROR);
    }
    l1file.spix = 0;
    l1file.epix = l1file.npix-1;
    l1file.dpix = 1;
    l1sensorID  = l1file.sensorID;
    
    if (l1sensorID == HMODISA || l1sensorID == HMODIST) {
        HY1C_out("-E- %s: The program does not currently support high resolution MODIS 16 bands.\n",argv[0]);
        exit(FATAL_ERROR);
    }
        
    
    nl1bands    = rdsensorinfo(l1sensorID, input.evalmask, NULL, NULL);
    i = rdsensorinfo(l1sensorID, input.evalmask, "Lambda", (void **) &l1wave);
    if ( (l1wavef = (float *)malloc(nl1bands*sizeof(float))) == NULL ) {
         HY1C_out("-E- %s: Error allocating memory to L1 wavelengths\n",argv[0]);
         exit(FATAL_ERROR);
    }
    for (i=0; i<nl1bands; i++) l1wavef[i] = (float)l1wave[i];
    nl1vbands   = windex(720.,l1wavef,(int)nl1bands);
    
    

    /*									*/
    /* Allocate memory for L1 scan data			  	        */
    /*									*/
    if ( alloc_l1(l1file.npix,NBANDS,NBANDSIR,l1file.n_inprods,&l1rec) == 0 ) {                                        
        HY1C_out("-E- %s: Unable to allocate L1 record.\n",argv[0]);
        exit(FATAL_ERROR);
    }
    
    /*									*/
    /* Allocate memory for L2/L3 data			      	        */
    /*									*/
    if (TRUE) {
    
    	sd_id = SDstart(input.il2file[0], DFACC_RDONLY);
        if (SDreadattr(sd_id,SDfindattr(sd_id,"Title"),(VOIDP)buf) == 0) {
	    
            if (SDreadattr(sd_id,SDfindattr(sd_id,"Sensor Name"),(VOIDP)buf) == 0) {
   
  		for (i=0; i<SENSOR_NUM; i++) {
       		    if (strcmp(buf, sensorName[i]) == 0 ) l2sensorID = i;
                }
	    }
	} else {
	     HY1C_out("-E- %s: Cannot recognize L3 HDF file\n",argv[0]);
	     exit(FATAL_ERROR);
	}
	// quick fix for TXC07
        if (l2sensorID == HMODISA) l2sensorID = MODISA;
        SDend(sd_id);
	
        nl2bands    = rdsensorinfo(l2sensorID, input.evalmask, NULL, NULL);
        i = rdsensorinfo(l2sensorID, input.evalmask, "Lambda", (void **) &l2wave);
	if ( (l2wavef = (float *)malloc(nl2bands*sizeof(float))) == NULL ) {
            HY1C_out("-E- %s: Error allocating memory to L2/L3 wavelengths\n",argv[0]);
            exit(FATAL_ERROR);
      	}
	for (i=0; i<nl2bands; i++) l2wavef[i] = (float)l2wave[i];
        nl2vbands   = windex(760.,l2wavef,(int)nl2bands);
	
	
	if (TRUE) {
  
	    if ( (fid = (int32 *)malloc(nl2files*sizeof(int32))) == NULL ) {
            	HY1C_out("-E- %s: Error allocating memory for reading L3 data.\n",argv[0]);
            	exit(FATAL_ERROR);
      	    }
	    if ( (sdfid = (int32 *)malloc(nl2files*sizeof(int32))) == NULL ) {
            	HY1C_out("-E- %s: Error allocating memory for reading L3 data.\n",argv[0]);
            	exit(FATAL_ERROR);
      	    }
 	    if ( (vgid = (int32 *)malloc(nl2files*sizeof(int32))) == NULL ) {
            	HY1C_out("-E- %s: Error allocating memory for reading L3 data.\n",argv[0]);
            	exit(FATAL_ERROR);
      	    }
 	    if ( (vdata_id = (l3nvdata *)malloc(nl2files*sizeof(l3nvdata))) == NULL ) {
            	HY1C_out("-E- %s: Error allocating memory for reading L3 data.\n",argv[0]);
            	exit(FATAL_ERROR);
      	    }
 	    if ( (meta_l3b = (meta_l3bType *)malloc(nl2files*sizeof(meta_l3bType))) == NULL ) {
            	HY1C_out("-E- %s: Error allocating memory for reading L3 data.\n",argv[0]);
            	exit(FATAL_ERROR);
      	    }
	    HY1C_out("%s: Matching L1 %s scene and L3 binned data from %s\n\n", argv[0], sensorName[l1sensorID], sensorName[l2sensorID]);

      	    for (i=0; i<nl2files; i++) {
	    
	    	strcpy(il3file, input.il2file[i]);
		status = open_input_l3b(il3file, &fid[i], &sdfid[i], &vgid[i], vdata_id[i], &meta_l3b[i], &n);
  	    	if (status == FAIL) {
        	    HY1C_out("-E- %s: Error opening L3 file %s.\n", argv[0], input.il2file[i]);
       		    exit(FATAL_ERROR);
		}
		if (i==0) nrows = n; else if (n != nrows) {
        	    HY1C_out("-E- %s: Mixing L3 files of different bin resolution, %d in %s and %d in %s\n", argv[0], nrows, input.il2file[0], n, input.il2file[i]);
       		    exit(FATAL_ERROR);
		}
		
  	    }
	}   
    } // end pixsrch==0

    /* set pixel and line limits */	
    if (input.epixl == -1)
        input.epixl = l1file.npix;
    if (input.eline == -1)
        input.eline = l1file.nscan;
    spix  = MAX(input.spixl - 1, 0);
    epix  = MIN(input.epixl - 1, l1file.npix-1);
    sscan = MAX(input.sline - 1, 0);
    escan = MIN(input.eline - 1, l1file.nscan-1);

    if (sscan > escan || spix > epix) {
        HY1C_out("-E- %s: scan and pixel limits make no sense.\n",argv[0]);
        HY1C_out(" start scan  = %d\n",sscan+1);
        HY1C_out(" end   scan  = %d\n",escan+1);
        HY1C_out(" start pixel = %d\n",spix +1);
        HY1C_out(" end   pixel = %d\n",epix +1);
        exit(FATAL_ERROR);
    }

    /* Note: for the L1 file, npix is still the native scan pixel count */
    l1file.spix = spix;                /* start pixel rel to L1 scan    */
    l1file.epix = epix;                /* end   pixel rel to L1 scan    */

    npix  = epix  - spix  + 1;
    nscan = escan - sscan + 1;


    if (TRUE) {
    
	if (TRUE) {
	

    	     /* Make sure number of L3 products are identical for every input L3 file */
    	     /* --------------------------------------------------------------------- */
     	     status = 0;
  	     nprod = MAXNVDATA-2;
  	     for (i=0; i<MAXNVDATA; i++) {
       		if (vdata_id[0][i] == -1) nprod--;
  	     }
    	     for (j=1; j<nl2files; j++) {
  	     	n = MAXNVDATA-2;
  	     	for (i=0; i<MAXNVDATA; i++) {
       		    if (vdata_id[j][i] == -1) n--;
  	        }
		if (n != nprod) {
		    HY1C_out("Number of products for %s (%d) differs from %s (%d)\n",
	            input.il2file[j], n, input.il2file[0], nprod);
		    status = -1;
    	     	}
    	     	if (status == -1) exit(FATAL_ERROR);
  	     }


    	     /* Make sure L3 product names are identical for every input L3 file */
    	     /* ---------------------------------------------------------------- */
    	     for (j=0; j<nl2files; j++) {

  	     	status = VSsetfields(vdata_id[j][0], "bin_num,nobs,weights,flags_set");
             	status = VSsetfields(vdata_id[j][1], "start_num,begin,extent,max");
	     }

  	     for (i=0; i<nprod; i++) {
	     
		VSgetname(vdata_id[0][2+i], data_prodname);
	        for (j=1; j<nl2files; j++) {
		
		    VSgetname(vdata_id[j][2+i], buf);
		    if (strcmp(buf, data_prodname) != 0) {
			HY1C_out("Product %d for %s (%s) differs from %s (%s)\n",
	        	i, input.il2file[j], buf, input.il2file[0], data_prodname);
			status = -1;
    	    	   }
    	           if (status == -1) exit(FATAL_ERROR);
		}
	     }


    	     /* Make sure all L2 products required for the cross-calibration are present */
    	     /* ------------------------------------------------------------------------ */ 
    	     chlinx = -1;
	     for (i=0; i<nprod; i++) {
		VSgetname(vdata_id[0][2+i], data_prodname);
       	   	if (strncmp("chlor_a\x0", data_prodname, 7) == 0) {
	  	     chlinx = i;
	      	     strcpy(buf, data_prodname);
	      	     strcat(buf, "_sum,");
	             strcat(buf, data_prodname);
	             strcat(buf, "_sum_sq");
		     for (j=0; j<nl2files; j++) status = VSsetfields(vdata_id[j][2+i], buf);
	  	     break;
       	        }
    	     }
    	     if (chlinx < 0) {
           	HY1C_out("-E- %s: Chlor_a product is missing from the L3 file data %s.\n",argv[0], input.il2file[0]);
           	exit(FATAL_ERROR);
    	     }
    	     nprods = 0;
    	     j = 0;
    	     for (i=0; i<nprod; i++) {
		VSgetname(vdata_id[0][2+i], data_prodname);
           	if (strncmp("nLw\x0", data_prodname, 3) == 0) {
	  	     prodinx[nprods++] = i;
	      	     strcpy(buf, data_prodname);
	      	     strcat(buf, "_sum,");
	             strcat(buf, data_prodname);
	             strcat(buf, "_sum_sq");
		     for (n=0; n<nl2files; n++) status = VSsetfields(vdata_id[n][2+i], buf);
	  	     j++;
           	}
	     }
    	     if (j < nl2vbands) {
	       HY1C_out("-E- %s: nLw products are missing from the L3 file data %s (%d of %d).\n", argv[0], input.il2file[0],j,nl2vbands);
           	exit(FATAL_ERROR);
    	     }
    	     if (j > nl2vbands) {
           	HY1C_out("-W- %s: nLw products for VIS and NIR bands in %s.\n", argv[0], input.il2file[0]);
		j = j - nl2vbands;
		nl1vbands += j;  
           	nl2vbands += j;
    	     }
    	     j = 0;
    	     aotinx = -1;
    	     for (i=0; i<nprod; i++) {
		VSgetname(vdata_id[0][2+i], data_prodname);
           	if (strncmp("tau\x0", data_prodname, 3) == 0) {
	  	     prodinx[nprods++] = i;
       	   	     if ((strncmp("taua_8\x0", data_prodname, 6) == 0) || (strncmp("tau_8\x0", data_prodname, 5) == 0)) aotinx = i;
	      	     strcpy(buf, data_prodname);
	      	     strcat(buf, "_sum,");
	             strcat(buf, data_prodname);
	             strcat(buf, "_sum_sq");
		     for (n=0; n<nl2files; n++) status = VSsetfields(vdata_id[n][2+i], buf);
	  	     j++;
       	   	}
	     }
    	     if (aotinx < 0) {
           	HY1C_out("-E- %s: AOT_865 product is missing from the L3 file data %s.\n", argv[0], input.il2file[0]);
           	exit(FATAL_ERROR);
    	     }
    	     if (j <= 1) {
           	HY1C_out("-W- %s: AOT products are missing from the L3 file data %s.\n", argv[0], input.il2file[0]);
           	nprods -= j;
		for (i=0; i<nl2bands; i++) prodinx[nprods++] = -1;
    	     }
    	     if (j > 1 && j < nl2bands) {
           	HY1C_out("-E- %s: AOT products are missing for some bands in %s.\n", argv[0], input.il2file[0]);
           	exit(FATAL_ERROR);
    	     }

  	     
	     for (i=0; i<nl2files; i++) {
	     
	     	VSinquire(vdata_id[i][2+prodinx[0]], &n_records, &interlace, fields, &vdata_size, vdata_name);
  
  		VSinquire(vdata_id[i][0], &n, &interlace0, fields, &vdata_size, vdata_name);
  		if (n != n_records) {
   		    HY1C_out("-E- %s: Different number of records for the product %s\n", argv[0], vdata_name);
        	    exit(FATAL_ERROR);
  		}
  		for (j=1; j<nprods; j++) {
       		    VSinquire(vdata_id[i][2+prodinx[j]], &n, &interlace, fields, &vdata_size, vdata_name);
       		    if (n != n_records) {
       	   		HY1C_out("-E- %s: Different number of records for the product %s\n", argv[0], vdata_name);
	   		exit(FATAL_ERROR);
       		    }
  		}
	     }
	     
	}
	
    } // end pixsrch==0
    
    
    minLon =  180.1;     /* min-longitude of L1 scene */
    maxLon = -180.1;   	 /* max-longitude of L1 scene */ 
    minLat =   90.1;   	 /* min-latitude  of L1 scene */
    maxLat =  -90.1;   	 /* max-latitude  of L1 scene */
    minLat1 =   90.1;    /* min-latitude  of L1 scene */
    maxLat1 =  -90.1;    /* max-latitude  of L1 scene */
    bndflg = 0;
    bndflg1 = 0;
    close2thepole = 0;
    close2thepole1 = 0;

    /*					 			*/
    /* Find min/max lat/lons for the L1 scene                   */
    /*								*/
    readl1(&l1file,sscan,&l1rec);
    mdate[0][0] = *(l1rec.year);
    mdate[0][1] = *(l1rec.day);
    mdate[0][2] = *(l1rec.msec);
    for (ipix=spix; ipix<=epix; ipix++) {
        if (minLat > l1rec.lat[ipix]) minLat = l1rec.lat[ipix];
    	if (maxLat < l1rec.lat[ipix]) maxLat = l1rec.lat[ipix];
    }
    if (minLat < -80.0) close2thepole = 1;
    if (maxLat > 80.0)  close2thepole = 1;
    if (close2thepole == 0) bndflg = (fabsf(l1rec.lon[spix] - l1rec.lon[epix]) > 180.0); else bndflg = -1;
    if (bndflg > 0) {
        /* first scan of L1 scene is crossing dateline */
        if (l1rec.lon[spix] >= 0.0 && minLon > l1rec.lon[spix]) minLon = l1rec.lon[spix]; 
	else if (l1rec.lon[epix] >= 0.0 && minLon > l1rec.lon[epix]) minLon = l1rec.lon[epix];
        if (l1rec.lon[epix] <= 0.0 && maxLon < l1rec.lon[epix]) maxLon = l1rec.lon[epix];
	else if (l1rec.lon[spix] <= 0.0 && maxLon < l1rec.lon[spix]) maxLon = l1rec.lon[spix];
    } else if (bndflg == 0) {
    	if (minLon > l1rec.lon[spix]) minLon = l1rec.lon[spix];
    	if (minLon > l1rec.lon[epix]) minLon = l1rec.lon[epix];
    	if (maxLon < l1rec.lon[spix]) maxLon = l1rec.lon[spix];
    	if (maxLon < l1rec.lon[epix]) maxLon = l1rec.lon[epix];
    }
    if (close2thepole == 0) {
    	if (bndflg == 0) 
	    if (l1rec.lon[spix] < l1rec.lon[epix]) west2east = 1; else west2east = 0;
	else
	    if (l1rec.lon[spix] > l1rec.lon[epix]) west2east = 1; else west2east = 0;
    }
    
    readl1(&l1file,escan,&l1rec);
    mdate[1][0] = *(l1rec.year);
    mdate[1][1] = *(l1rec.day);
    mdate[1][2] = *(l1rec.msec);
    for (ipix=spix; ipix<=epix; ipix++) {
        if (minLat1 > l1rec.lat[ipix]) minLat1 = l1rec.lat[ipix];
    	if (maxLat1 < l1rec.lat[ipix]) maxLat1 = l1rec.lat[ipix];
    }
    if (minLat1 < -80.0) close2thepole1 = 1;
    if (maxLat1 > 80.0)  close2thepole1 = 1;
    if (close2thepole1 == 0) bndflg1 = (fabsf(l1rec.lon[spix] - l1rec.lon[epix]) > 180.0); else bndflg1 = -1;
    if (bndflg1 > 0) {
        /* last scan of L1 scene is crossing dateline */
        if (l1rec.lon[spix] >= 0.0 && minLon > l1rec.lon[spix]) minLon = l1rec.lon[spix]; 
	else if (l1rec.lon[epix] >= 0.0 && minLon > l1rec.lon[epix]) minLon = l1rec.lon[epix];
        if (l1rec.lon[epix] <= 0.0 && maxLon < l1rec.lon[epix]) maxLon = l1rec.lon[epix];
	else if (l1rec.lon[spix] <= 0.0 && maxLon < l1rec.lon[spix]) maxLon = l1rec.lon[spix];
    } else if (bndflg1 == 0) {
    	if (minLon > l1rec.lon[spix]) minLon = l1rec.lon[spix];
    	if (minLon > l1rec.lon[epix]) minLon = l1rec.lon[epix];
    	if (maxLon < l1rec.lon[spix]) maxLon = l1rec.lon[spix];
    	if (maxLon < l1rec.lon[epix]) maxLon = l1rec.lon[epix];
    }
    if (close2thepole1 == 0) {
    	if (bndflg1 == 0) 
	    if (l1rec.lon[spix] < l1rec.lon[epix]) west2east1 = 1; else west2east1 = 0;
	else
	    if (l1rec.lon[spix] > l1rec.lon[epix]) west2east1 = 1; else west2east1 = 0;
    }
    
    if (close2thepole == 0 && close2thepole1 == 0) {
    	if (bndflg == 0 && bndflg1 == 0) {
	    if ((west2east && west2east1 == 0) || (west2east1 && west2east == 0)) {
	    	close2thepole = 1;
		close2thepole1 = 1;
	    }
	} else if ((bndflg && bndflg1 == 0) || (bndflg1 && bndflg == 0)) {
	    if ((west2east && west2east1) || (west2east1 == 0 && west2east == 0)) {
	    	close2thepole = 1;
		close2thepole1 = 1;
	    }
	}
    }
    	
    if (close2thepole == 0 && close2thepole1 == 0) {
        if (minLat1 < minLat) minLat = minLat1;
        if (maxLat1 > maxLat) maxLat = maxLat1;
    } else
    if (close2thepole1 == 0 && close2thepole) {
        if (minLat < 0.0) {
	    minLat = -90.0;
	    maxLat = maxLat1;
	} else 
	if (maxLat > 0.0) {
	    maxLat = 90.0;
	    minLat = minLat1;
	}
    } else
    if (close2thepole == 0 && close2thepole1) {
        if (minLat1 < 0.0) minLat = -90.0; else maxLat = 90.0;
    } else
    if (close2thepole && close2thepole1) {
        if (minLat1 < 0.0) minLat = -90.0; else if (maxLat1 > 0.0) maxLat = 90.0;
        if (minLat < 0.0) minLat = -90.0; else if (maxLat > 0.0) maxLat = 90.0;
    }
    
    if (bndflg > 0 || bndflg1 > 0) bndflg = 1; else bndflg = 0;
    if (close2thepole || close2thepole1) close2thepole = 1;
    meanday(mdate[0][0], mdate[0][1], mdate[0][2], mdate[1][0], mdate[1][1], mdate[1][2], &(ldate[0]), &(ldate[1]), &(ldate[2]));
    

    HY1C_out("\nL1 data:\n");
    HY1C_out("Min Lat = %f, Max Lat = %f, Min Lon = %f, Max Lon = %f\n", minLat, maxLat, minLon, maxLon);
    HY1C_out("spixl = %d, epixl = %d, sscan = %d, escan = %d\n", spix, epix, sscan, escan);
    HY1C_out("year = %d, day = %d, msec = %d\n\n", ldate[0], ldate[1], ldate[2]);

        
    
    /*					 			*/
    /* Check which are the L2/L3 files that overlap the L1 scene   */
    /*					 			*/
    if (TRUE) {

    	if ( (l2finx = (int *)malloc(nl2files*sizeof(int))) == NULL ) {
           HY1C_out("-E- %s: Error allocating memory for L2 file index.\n",argv[0]);
           exit(FATAL_ERROR);
    	}
    
	if (TRUE) {
    
           if ( (bins = (int32 **)malloc(nl2files*sizeof(int32 *))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory to bin numbers\n",argv[0]);
           	exit(FATAL_ERROR);
           }
	   
           if ( (nbins = (int32 *)malloc(nl2files*sizeof(int32))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory to number of bins\n",argv[0]);
           	exit(FATAL_ERROR);
           }

    	   if ( (flagusemask = (uint32_t *)malloc(1*sizeof(uint32_t))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory for L2 flags.\n",argv[0]);
           	exit(FATAL_ERROR);
    	   }
	
           if ( (sum_prod = (float32 ***) calloc(nl2files, sizeof(float32 **))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory to L3 file products\n",argv[0]);
           	exit(FATAL_ERROR);
	   }
           if ( (sum_chl  = (float32 **) calloc(nl2files, sizeof(float32 *))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory to L3 file chlor_a\n",argv[0]);
           	exit(FATAL_ERROR);
	   }

  	   /* ----------------- */
  	   /* Set up bin arrays */
  	   /* ----------------- */
  	   if ( (numbin   = (int32 *)calloc(nrows, sizeof(int32))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory to numbin\n",argv[0]);
           	exit(FATAL_ERROR);
  	   }
  	   if ( (basebin  = (int32 *)calloc(nrows+1, sizeof(int32))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory to basebin\n",argv[0]);
           	exit(FATAL_ERROR);
  	   }

  	   for (i=0; i<nrows; i++) {
    	   	f32 = (i + 0.5) * (180.0 / nrows) - 90.0;
     	   	numbin[i] = (int32) (cos(f32 * PI/ 180.0) * (2.0*nrows) + 0.5);
           }

  	   basebin[0] = 1;
  	   for (i=1; i<=nrows; i++) {
    	   	basebin[i] = basebin[i-1] + numbin[i-1];
  	   }
	
	}
	

    	n = 0;
        for (i=0; i<nl2files; i++)  {

    
	   l2finx[i] = 1;

	   if (TRUE) {
	   
           	/* grab the target sensor name from the L3 meta-data */
	   	if ( i==0 ) 
                    strcpy(l2sensor, meta_l3b[i].sensor_name); 
           	else if (strcmp(l2sensor, meta_l3b[i].sensor_name) != 0) {
                    HY1C_out("-E- %s: Mixing L3 data from different sensors, %s and %s.\n",
                 	argv[0], l2sensor, meta_l3b[i].sensor_name);
	            exit(FATAL_ERROR);
    	   	}
	   
	   
	   	for (j=0; j<i; j++) {
	        if (meta_l3b[i].syear == meta_l3b[j].syear && meta_l3b[i].sday ==
	     	   meta_l3b[j].sday && meta_l3b[i].smsec == meta_l3b[j].smsec &&
	     	   meta_l3b[i].eyear == meta_l3b[j].eyear && meta_l3b[i].eday ==
	     	   meta_l3b[j].eday && meta_l3b[i].emsec == meta_l3b[j].emsec) {
	     	  	HY1C_out("-W- %s: L2 coverage is repeated between the files %s and %s\n\n", 
		  	argv[0], input.il2file[i], input.il2file[j]);
		  	l2finx[i] = 0;
	           }
	   	}

                if ((yds2unix(meta_l3b[i].eyear, meta_l3b[i].eday, meta_l3b[j].emsec/1000) -
                     yds2unix(meta_l3b[i].syear, meta_l3b[i].sday, meta_l3b[j].smsec/1000)) > (time_t) 86400 * 365) {
		    climatology=1;
                    HY1C_out("Assuming climatology, disabling coefficient of variation test.\n");
		}

	
           	/* Setup flag mask */
           	/* --------------- */
/*           	strcpy(buf, meta_l3b[i].flag_names);   */
           	strcpy(buf, "\x0");
           	for (j=0; j<NFLAGS; j++) {
               	   strcat(buf, l2_flag_lname[j]);
               	   if (j < NFLAGS-1) strcat(buf, ",\x0");
          	}
           	setupflags(buf, input.flaguse, &(flagusemask[0]), &required_mask, &status );
           	if (status < 0) {
         	    HY1C_out("-E- %s: Error reading L2 flags %s.\n",argv[0], input.il2file[i]);
         	    exit(FATAL_ERROR);
           	}

	   
	   	VSinquire(vdata_id[i][0], &n_records, &interlace0, fields, &vdata_size, vdata_name);
	   	VSinquire(vdata_id[i][2+prodinx[0]], &n_records, &interlace, fields, &vdata_size, vdata_name);
	   
           	if ( (bins[i] = (int32 *)malloc(n_records*sizeof(int32))) == NULL ) {
         	    HY1C_out("-E- %s: Error allocating memory to bin numbers\n",argv[0]);
         	    exit(FATAL_ERROR);
           	}
	   	   
  	   	vsize = 2*sizeof(int32)+sizeof(float32)+sizeof(int16);
  	   	if ( (binlist_buf = (uint8 *) calloc(n_records, vsize)) == NULL ) {
        	    HY1C_out("-E- %s: Error allocating memory to L3 file bin info\n",argv[0]);
        	    exit(FATAL_ERROR);
  	   	}
  	   	for (j=0; j<2; j++) {
        	    if ( (sum_buf[j] = (float32 *) calloc(2*n_records, sizeof(float32))) == NULL ) {
            	    	HY1C_out("-E- %s: Error allocating memory to L3 file products\n",argv[0]);
            	    	exit(FATAL_ERROR);
		    }
  	   	} /* iprod loop */
        	if ( (sum_chl[i] = (float32 *) calloc(2*n_records, sizeof(float32))) == NULL ) {
            	    HY1C_out("-E- %s: Error allocating memory to L3 file chlor_a\n",argv[0]);
            	    exit(FATAL_ERROR);
		}
        	if ( (sum_prod[i] = (float32 **) calloc(nprods, sizeof(float32 *))) == NULL ) {
            	    HY1C_out("-E- %s: Error allocating memory to L3 file products\n",argv[0]);
            	    exit(FATAL_ERROR);
		}
  	   	for (j=0; j<nprods; j++) {
        	    if ( (sum_prod[i][j] = (float32 *) calloc(2*n_records, sizeof(float32))) == NULL ) {
            	    	HY1C_out("-E- %s: Error allocating memory to L3 file products\n",argv[0]);
            	    	exit(FATAL_ERROR);
		    }
  	   	} /* iprod loop */
			   

      		VSseek(vdata_id[i][0], 0);
  	   	binindex_buf = n_records;
  	   	/* Get bin numbers for each filled bin */
 	   	if (VSread(vdata_id[i][0], (uint8 *) binlist_buf, binindex_buf, interlace0) == -1) {
      		    HY1C_out("-E- %s: Unable to read bin numbers from %s\n", argv[0], input.il2file[i]);
      		    exit(FATAL_ERROR);
  	   	}
  	   	/* error check */
  	   	if (binindex_buf != n_records) {
      		    HY1C_out("-E- %s: Wrong number of records were read: %d\n", argv[0], binindex_buf);
      		    exit(FATAL_ERROR);
  	   	}
	   
	   	if (input.chlthreshold < CHL_MAX || input.aotthreshold < AOT_MAX) {
		
      		    if (input.chlthreshold < CHL_MAX) {
      			VSseek(vdata_id[i][2+chlinx], 0);
		    	if (VSread(vdata_id[i][2+chlinx], (uint8 *) sum_buf[0], binindex_buf, interlace) == -1) {
      			    HY1C_out("-E- %s: Unable to read chlorophyll binned data from %s\n", argv[0], input.il2file[i]);
      			    exit(FATAL_ERROR);
  	   	    	}
  	   	    	/* error check */
  	   	    	if (binindex_buf != n_records) {
      			    HY1C_out("-E- %s: Wrong number of records were read: %d\n", argv[0], binindex_buf);
      			    exit(FATAL_ERROR);
  	   	    	}
	   	    }
      		    if (input.aotthreshold < AOT_MAX) {
      			VSseek(vdata_id[i][2+aotinx], 0);
		        if (VSread(vdata_id[i][2+aotinx], (uint8 *) sum_buf[1], binindex_buf, interlace) == -1) {
      			    HY1C_out("-E- %s: Unable to read AOT binned data from  %s\n", argv[0], input.il2file[i]);
      			    exit(FATAL_ERROR);
  	   	    	}
  	   	    	/* error check */
  	   	    	if (binindex_buf != n_records) {
      			    HY1C_out("-E- %s: Wrong number of records were read: %d\n", argv[0], binindex_buf);
      			    exit(FATAL_ERROR);
  	   	    	}
	   	    }
	   	}
		
  		/* Get data values (sum) for each filled bin */
  		/* ----------------------------------------- */
  		for (j=0; j<nprods; j++) {
  

      		    VSseek(vdata_id[i][2+prodinx[j]], 0);
      
      		    if (VSread(vdata_id[i][2+prodinx[j]], (uint8 *) sum_prod[i][j], binindex_buf, interlace) == -1) {
	  		HY1C_out("-E- %s: Unable to read sum/sum_sq: %d\n", argv[0], binindex_buf);
          		exit(FATAL_ERROR);
      		    }
      
       		    /* error check */
      		    if (binindex_buf != n_records) {
      	  		HY1C_out("-E- %s: Wrong number of records were read: %d\n", argv[0], binindex_buf);
          		exit(FATAL_ERROR);
      		    }
      
  		} /* iprod loop */
  
	   
  		close_input_l3b(fid[i], sdfid[i], vgid[i], vdata_id[i]);
		
		
                if( input.vct_msk_water == 1 )
                  {
        	  if ( (mask = (char *) malloc(5940424 * sizeof(char))) == NULL ) {
            	      HY1C_out("-E- %s: Error allocating memory to the mask data\n", argv[0]);
            	      exit(FATAL_ERROR);
		    }
  		  strcpy(maskname, input.vct_water );
  	
  		  if (read9km_mask(maskname, mask) > 0) {
              	      HY1C_out("-E- %s: Failed reading the deep water mask file\n", argv[0]);
              	      exit(FATAL_ERROR);
          	  }
		}
		   	   
  	   	no = 0;
		wgtsum = 0.0;
  	   	for (j=0; j<n_records; j++) {
      
      		    k = 1;
		    memcpy((void *)&bin_num,(const void *)&binlist_buf[j*vsize],    4);
      		    memcpy((void *)&nobs,   (const void *)&binlist_buf[j*vsize+4],  2);
      		    memcpy((void *)&wgt,    (const void *)&binlist_buf[j*vsize+6],  4);
      		    memcpy((void *)&flags,  (const void *)&binlist_buf[j*vsize+10], 4);
		    
                    if( ( input.vct_msk_water == 1 ) && 
                      ( (int32_t)is_masked(bin_num, mask, nrows) <= 0 ) ) 
                      continue;
		    
      		    if (nobs <= 0 || wgt <= 0.0 || 
                      (flags & flagusemask[0])) continue;
		   /* check for min # samples in a bin for data use */
                    if( nobs < input.vct_min_nbin ) continue;
		    
		    if (input.aotthreshold < AOT_MAX) {
		     	memcpy((void *)&f32, (const void *)&sum_buf[1][2*j], 4);
	                memcpy((void *)&v32, (const void *)&sum_buf[1][2*j+1], 4);
			f32 /= wgt;
  	      	        v32 = sqrt(v32/wgt - f32*f32);
		     	if (f32 > input.aotthreshold) k = 0;
                        if (!climatology)
			  if (v32 / f32 > 0.15) k = 0;
		    }
		    if (input.chlthreshold < CHL_MAX) {
		     	memcpy((void *)&f32, (const void *)&sum_buf[0][2*j], 4);
	                memcpy((void *)&v32, (const void *)&sum_buf[0][2*j+1], 4);
			f32 /= wgt;
  	      	        v32 = sqrt(v32/wgt - f32*f32);
		     	if (f32 > input.chlthreshold) k = 0;
                        if (!climatology)
  			  if (v32 / f32 > 0.15) k = 0;
		    }
		    
      		    if (k <= 0) continue;
		    
		    
		    l = -1;
		    while (bin_num >= basebin[l+1]) l++;
      		    lat = ((l + 0.5) / nrows) * 180.0 - 90.0;
      		    lon = 360.0 * (bin_num - basebin[l] + 0.5) / numbin[l] - 180.0;
		   
		
	  	    bin_row = (int32_t) ((90.0 + lat) * ((float64) nrows / (float64) 180.0));
	  	    bin_col = (int32_t) ((float64) numbin[bin_row] * (lon + 180.0) / (float64) 360.0);
	  	    bin = basebin[bin_row] + bin_col;
		    
		    if (bin_num != bin) {
	  		HY1C_out("-E- %s: Lat/lon wrongly calculated %f, %f\n", argv[0], lat, lon);
          		exit(FATAL_ERROR);
      		    }
		    
		    /*  check which L3 bins lie within the L1 coverage  */
		    latTest = (lat <= maxLat+0.2 && lat >= minLat-0.2);
	   	    if (close2thepole) {
		     	lonTest = 1;
	   	    }
		    else {
		  	if (bndflg > 0) {
	    	    	    lonTest = (lon > maxLon+0.2 && lon < minLon-0.2);
		    	    lonTest = 1 - lonTest; 
		     	} else 
	    	    	    lonTest = (lon <= maxLon+0.2 && lon >= minLon-0.2);
		    }
		    if (lonTest + latTest != 2) k = 0;
		    
      		    if (k <= 0) continue;
		    	    
		    
      		    if (nobs > 0 && wgt > 0.0 && ((flags & flagusemask[0]) == 0) && k > 0) {
	  	     	bins[i][no] = bin_num;
			wgtsum += wgt;
  			for (l=0; l<nprods; l++) sum_prod[i][l][no] = sum_prod[i][l][j*2]/wgt;
			sum_chl[i][no] = f32;
			no++;
		    }			
  	        }
	   
	   	nbins[i] = no;
  	   	free(binlist_buf);
  	   	for (j=0; j<2; j++) free(sum_buf[j]);
  	   	if( input.vct_msk_water == 1 ) 
                  {
                  free(mask);
	   	  j = (int32_t)is_masked(bin_num, mask, -1);
                  }
	   
  		if (no <= 0) {
		
		    HY1C_out("-E- %s: There are no valid bins in file %s\n",argv[0], input.il2file[i]);
		    free(bins[i]);
		    for (j=0; j<nprods; j++) free(sum_prod[i][j]);
		    free(sum_prod[i]);
		    free(sum_chl [i]);
		    l2finx[i] = 0;
		    
		} else {
		
		    HY1C_out("Mean weights per file %d, %f\n", i+1, wgtsum/nbins[i]);
		    
  	   	    if (no < n_records) {
		
       		        if ((realloc((void *)bins[i], nbins[i]*sizeof(int32))) == NULL ) {
            	    	    HY1C_out("-E- %s: Error reallocating memory to L3 file bins\n",argv[0]);
            	    	    exit(FATAL_ERROR);
       		        }
       		        for (j=0; j<nprods; j++) if ((realloc((void *)sum_prod[i][j], nbins[i]*sizeof(float32))) == NULL ) {
            	    	    HY1C_out("-E- %s: Error reallocating memory to L3 file products\n",argv[0]);
            	    	    exit(FATAL_ERROR);
       		        }
       		        if ((realloc((void *)sum_chl[i], nbins[i]*sizeof(float32))) == NULL ) {
            	    	    HY1C_out("-E- %s: Error reallocating memory to L3 file chlor_a\n",argv[0]);
            	    	    exit(FATAL_ERROR);
       		        }
  	   	    }
		}
		
                HY1C_out("Number of valid data bins in the %d L3 file: %d\n", i+1, nbins[i]);

	   }
	   
           n += l2finx[i];

       	}


	   
        free(flagusemask);
	free(fid);
	free(sdfid);
	free(vgid);
 	free(vdata_id); 

	
    	/*					 			    */
        /* if no valid L2/L3 files to target, then bail                        */
    	/*					 			    */
    	if ( n < 1 ) {
	
	    free(meta_l3b); 
            free(bins);
	    free(nbins);
	    free(sum_prod);
	    free(sum_chl);
	    free(numbin);
	    free(basebin);
	    
    	    closel1(&l1file);
    	    free_l1(&l1rec);
	    free(l2finx);
	    free(l1wavef);
	    free(l2wavef);

           /*
            *  make an empty target, aer file set if desired
            */
            if( input.vct_fill_pct < 0. )
              empty_tgtaer( npix, nscan, l2sensorID, outtarget, outaer );

            HY1C_out("-E- %s: No valid L2/L3 ocean data within the L1 scene %s\n", argv[0], input.ifile[0]);
	    
            exit(NOMATCH_ERROR);
    	}
	

    
	
    	/*					 			    */
    	/* Free some memory if not all L2 files overlap the L1 scene        */
    	/*					 			    */
    	if ( n < nl2files ) {
	
    	    if ( (inx = (int32_t *)malloc(n*sizeof(int32_t))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory for L2 file index.\n",argv[0]);
           	exit(FATAL_ERROR);
    	    }
	    j = 0;
	    for (i=0; i<nl2files; i++) 
	    	if (l2finx[i]) {
	            inx[j++] = i;
	        }
		
	    nl2files = n;
    	    if ( (l2finx = (int *)realloc(l2finx, nl2files*sizeof(int))) == NULL ) {
             	HY1C_out("-E- %s: Error reallocating memory for L2 file index.\n",argv[0]);
             	exit(FATAL_ERROR);
            }
       	}
	else {
	
    	    if ( (inx = (int32_t *)malloc(nl2files*sizeof(int32_t))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory for L2 file index.\n",argv[0]);
           	exit(FATAL_ERROR);
    	    }
	    for (i=0; i<nl2files; i++) inx[i] = i;
	}
	
	
    	/*					 			    */
    	/* Order L2 files depending on their distance in time from L1 scene */
    	/*					 			    */
	/*
	if (nl2files > 1) {
	
    	    if ( (timediff = (uint64_t *)malloc(nl2files*sizeof(uint64_t))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory for time difference vector.\n",argv[0]);
           	exit(FATAL_ERROR);
    	    }
    	    if ( (xtimediff = (uint64_t *)malloc(nl2files*sizeof(uint64_t))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory for time difference vector.\n",argv[0]);
           	exit(FATAL_ERROR);
    	    }
    	    if ( (d2s = (int32_t *)malloc(nl2files*sizeof(int32_t))) == NULL ) {
           	HY1C_out("-E- %s: Error allocating memory for L2 file pixel distances.\n",argv[0]);
           	exit(FATAL_ERROR);
    	    }
	    for (i=0; i<nl2files; i++) {
	        meanday((int32_t)meta_l3b[inx[i]].syear, (int32_t)meta_l3b[inx[i]].sday, meta_l3b[inx[i]].smsec, (int32_t)meta_l3b[inx[i]].eyear, (int32_t)meta_l3b[inx[i]].eday, meta_l3b[inx[i]].emsec, &(mdate[0][0]), &(mdate[0][1]), &(mdate[0][2]));
	   	timediff[i] = diffmsec(ldate[0], ldate[1], ldate[2], mdate[0][0], mdate[0][1], mdate[0][2]);
		xtimediff[i] = timediff[i];
	    }
	    
	    qsort((void *)xtimediff, nl2files, sizeof(uint64_t), &compar);

	    status = 0;
	    for (i=0; i<nl2files; i++) {
	    	for (j=0; j<nl2files; j++) if (xtimediff[i] == timediff[j]) d2s[i] = inx[j];
	    	if (d2s[i] != inx[i]) status = 1;
	    }
	    
	    if ( status ) {
	    
	    	for (i=0; i<nl2files; i++) inx[i] = d2s[i]; 
	    }
	     
	    free(timediff);
	    free(xtimediff);
	    free(d2s);
	     
    	    HY1C_out("There are %d L2/L3 files overlapping the L1 scene\n", nl2files);
		    
       	}
	else
    	    HY1C_out("There is %d L2/L3 file overlapping the L1 scene\n", nl2files);
	*/


    } // end pixsrch==0



    /*					 			    */
    /* Initialize scan line output variables                        */
    /*					 			    */
    
    if ((tglength = alloc_target(npix,NBANDS,&tgtstr)) == 0) {
        HY1C_out("%s: Error: Unable to allocate memory for L2 binary target file data %s.\n", argv[0], outtarget);
        exit(FATAL_ERROR);
    }
    if ((arlength = alloc_aer(npix,NBANDS,&aerstr)) == 0) {
        HY1C_out("%s: Error: Unable to allocate memory for L2 binary aerosol file data %s.\n", argv[0], outaer);
        exit(FATAL_ERROR);
    }

    /*					 			    */
    /* Initialize the header for the L2 binary output               */
    /*					 			    */
    tgtstr.sensorID = l2sensorID;
    l2hdr[0] = l2sensorID;
    l2hdr[1] = tglength;
    l2hdr[2] = npix;
    l2hdr[3] = 0;
    l2hdr[4] = nscan;
    
    
    /*					 			    */
    /* Allocate space for target and aerosol record blocks          */
    /*					 			    */
    tout = (unsigned char *)malloc(OUTSCANLINE*tglength);
    aout = (unsigned char *)malloc(OUTSCANLINE*arlength);
    if (tout == NULL || aout == NULL) { 
         HY1C_out("%s: Error: Unable to allocate memory for the %d scan-line L2/L1 binary output %s %s, exiting...\n", argv[0], OUTSCANLINE, outtarget, outaer);
         exit(FATAL_ERROR);
    }
    scan100 = OUTSCANLINE;
    pcvalid = 0L;
    badconv = 0L;
    ff  = nl1bands*sizeof(float);
    ll = 3*sizeof(int32_t);
    
    initialTime = time(&initialTime);


    if (TRUE) {

    
    	if ( (nodata = (float *)malloc(NBANDS*sizeof(float))) == NULL ) {
            HY1C_out("-E- %s: Error allocating memory to the zeroed vector.\n",argv[0]);
            exit(FATAL_ERROR);
    	}
        for (j=0; j<NBANDS; j++) nodata[j] = -1.0;

    	/*					 		      */
    	/* Read L1 file scan by scan, find L2 data and write binaries.*/
    	/*							      */
    	oscan = 0;
	wscan = 0;
        interp_swf_mod = 0;
        if( ( l1sensorID == MODISA || l1sensorID == MODIST) && 
          ( l2sensorID == SEAWIFS ) ) interp_swf_mod = 1;

    	for (iscan = sscan; iscan <= escan; iscan++) {


      	    if (readl1(&l1file,iscan,&l1rec) != 0) {
            	fHY1C_out(stderr,
            	"-E- %s Line %d: error reading %s at scan %d.\n",
            	__FILE__,__LINE__,l1file.name,iscan);
            	exit(FATAL_ERROR);
      	    }

      	    if (iscan % 50 == 0) {
                finalTime = time(&finalTime);
           	HY1C_out("scan %d (lat %f), out of %d scans, time from start %f [min]\n", iscan, l1rec.lat[npix/2], escan+1, (float)(finalTime-initialTime)/60.0);
	    }
	    
      	    for (i=0; i<3; i++) udate[i] = 0;
	    
      	    fpixs = 0;
      	    opix = 0;
	     

      	    for (ipix = spix; ipix <= epix; ipix++) {
	

	  	lat = l1rec.lat[ipix];
	  	lon = l1rec.lon[ipix];
	  
	  	if (lon > 180) lon = lon - 360;
	  	status = 0;
		
	  
		n = 0;


	  	if (TRUE) {
		
	  	   bin_row = (int32_t) ((90.0 + lat) * ((float64) nrows / (float64) 180.0));
	  	   bin_col = (int32_t) ((float64) numbin[bin_row] * (lon + 180.0) / (float64) 360.0);
	  	   bin = basebin[bin_row] + bin_col;

		   for (i=0; i<nl2files; i++) {

		       ibin = binlocate(bins[inx[i]], nbins[inx[i]], bin);

		       if (ibin > 0) {
		    
		            l = -1;
		            while (bins[inx[i]][ibin] >= basebin[l+1]) l++;
      		            lat1 = ((l + 0.5) / nrows) * 180.0 - 90.0;
      		            lon1 = 360.0 * (bins[inx[i]][ibin] - basebin[l] + 0.5) / numbin[l] - 180.0;

		            wgtsum = sqrt((lat-lat1)*(lat-lat1)+(lon-lon1)*(lon-lon1));
		            if (wgtsum > 0.1) {
         		    	HY1C_out("%s: Warning: The difference in lat/lons is large, %f.\n", argv[0], wgtsum);
		            }

		  	    chl = sum_chl[inx[i]][ibin];
			    
			    tgtstr.solz[opix] = -1.0;
			    
			    for (j=0; j<nl1vbands; j++) {
                                if( interp_swf_mod == 0 )
                                  ib = ( j < nl2vbands ) ? j : nl2vbands - 1;
                                else
                                  ib = windex(l1wavef[j],l2wavef,(int)nl2bands);
			    	tgtstr.Lw[opix*NBANDS+j] = -1.0;
			    	if (prodinx[ib] >= 0) tgtstr.nLw[opix*NBANDS+j] = sum_prod[inx[i]][ib][ibin]; else tgtstr.nLw[opix*NBANDS+j] = -1.0;
			    }
			    for (j=nl1vbands; j<NBANDS; j++) {
			    	tgtstr.Lw[opix*NBANDS+j] = -1.0;
			    	if (j < nl1bands) tgtstr.nLw[opix*NBANDS+j] = 0.0; else tgtstr.nLw[opix*NBANDS+j] = -1.0;
			    }
			    
	  		    aerstr.mod_min[opix] = -1;
	  		    aerstr.mod_max[opix] = -1;
	  		    aerstr.mod_rat[opix] = -1.0;
	  		    
			    for (j=0; j<nl1bands; j++) {
                                if( interp_swf_mod == 0 )
                                  ib = ( j < nl2vbands ) ? j : nl2vbands - 1;
                                else
                                  ib = windex(l1wavef[j],l2wavef,(int)nl2bands);
			    	if (prodinx[nl2vbands+ib] >= 0) aerstr.taua[opix*NBANDS+j] = sum_prod[inx[i]][nl2vbands+ib][ibin]; else aerstr.taua[opix*NBANDS+j] = -1.0;
			    }				
			    for (j=nl1bands; j<NBANDS; j++) {
			    	aerstr.taua[opix*NBANDS+j] = -1.0;
			    }				
			    
      		  	    udate[0] += (uint64_t)meta_l3b[inx[i]].syear;
      		  	    udate[1] += (uint64_t)meta_l3b[inx[i]].sday;
      		  	    udate[2] += (uint64_t)meta_l3b[inx[i]].smsec;
       		  	    udate[0] += (uint64_t)meta_l3b[inx[i]].eyear;
       		  	    udate[1] += (uint64_t)meta_l3b[inx[i]].eday;
       		  	    udate[2] += (uint64_t)meta_l3b[inx[i]].emsec;
 		  	    ++fpixs;
			    ++n;
			    
     			    break;
			    
			}
			   
		   }

	  	}

		
	  	if ( n > 0 ) {
	          /*  Note that for all other cases than SeaWiFS2MODIS,
                      a spline interpolation is applied to all available 
                      bands, else, use a chl based interp to nearest bands 
                   */			
                   if( interp_swf_mod == 1 ) {
		       status = nLw_SeaWiFS2MODIS(l1sensorID, &(tgtstr.nLw[opix*NBANDS]), &(tgtstr.Lw[opix*NBANDS]));
		       if (status == 1) {
		           ++badconv;
			   n = 0L;
		       }
		   } else {
		       bandpass_correct(l2sensorID,l1sensorID,l2wavef,l1wavef,nl2bands,nl1bands,&(tgtstr.nLw[opix*NBANDS]),&(tgtstr.Lw[opix*NBANDS]));
		   }
                   if (l1sensorID != l2sensorID)
                       tau_interp(nl1bands, l1wavef, nl2bands, l2wavef, &(aerstr.taua[opix*NBANDS]));
		   if ( n > 0 ) ++pcvalid;
	  	}
		
		if ( n < 1 ) {
		   tgtstr.solz[opix] = -1.0;
		   memcpy((void *) &(tgtstr.Lw[opix*NBANDS]),  (const void *) nodata, ff); 
		   memcpy((void *) &(tgtstr.nLw[opix*NBANDS]), (const void *) nodata, ff);
		   aerstr.mod_min[opix] = -1L;
	  	   aerstr.mod_max[opix] = -1L;
	  	   aerstr.mod_rat[opix] = -1.0;
		   memcpy((void *) &(aerstr.taua[opix*NBANDS]), (const void *) nodata, ff);
	  	}
		
	  	++opix;
	   
      	    }
      
      
            if (fpixs > 0) {
      		for (j=0; j<3; j++) {
		   udate[j] /= (uint64_t)(2*fpixs);
		   ldate[j] = (int32_t)udate[j];
		}
	        memcpy((void *) tgtstr.data, (const void *) ldate, ll);
            }
	    else {
	        *tgtstr.year = 0L;
		*tgtstr.day  = 0L;
		*tgtstr.msec = 0L;
	    }
	    
      	    if (opix != npix) {
        	HY1C_out("%s: Error: Incorrect number of pixels %d obtained in the scan %d.\n", argv[0], opix, oscan);
            }
	    
	    memcpy((void *) &(tout[wscan*tglength]), (const void *)(tgtstr.data), tglength);
	    memcpy((void *) &(aout[wscan*arlength]), (const void *)(aerstr.data), arlength);

            ++oscan;
	    ++wscan;
	    
	    
    	    /* Write out if data sent in blocks of OUTSCANLINE scans   */
	    
	    if (oscan % scan100 == 0) {
	    
    		if (oscan == OUTSCANLINE) {
		    if ((outtgt = fopen(outtarget,"w")) == NULL) {
        	    	HY1C_out("%s: Error: Unable to open %s for writing.\n", argv[0], outtarget);
        	    	exit(FATAL_ERROR);
    		    }
    		    l2hdr[1] = tglength;
    		    fwrite(l2hdr,sizeof(l2hdr),1,outtgt);
    		    fseek(outtgt,tglength,SEEK_SET);
		} else
		    if ((outtgt = fopen(outtarget,"a")) == NULL) {
        	    	HY1C_out("%s: Error: Unable to open %s for appending.\n", argv[0], outtarget);
        	    	exit(FATAL_ERROR);
    		    }
    		fwrite((const void *)tout, OUTSCANLINE*tglength, 1, outtgt);
    		fclose(outtgt);
	
    		if (oscan == OUTSCANLINE) {
    		    if ((outar = fopen(outaer,"w")) == NULL) {
        		HY1C_out("%s: Error: Unable to open %s for writing.\n", argv[0], outaer);
        		exit(FATAL_ERROR);
    		    }
    		    l2hdr[1] = arlength;
    		    fwrite(l2hdr,sizeof(l2hdr),1,outar);
    		    fseek(outar,arlength,SEEK_SET);
		} else
     		    if ((outar = fopen(outaer,"a")) == NULL) {
        		 HY1C_out("%s: Error: Unable to open %s for appending.\n", argv[0], outaer);
        		 exit(FATAL_ERROR);
		    }
   		fwrite((const void *)aout, OUTSCANLINE*arlength, 1, outar);
    		fclose(outar);
		
		wscan = 0L;
	    }
	    
 
    	}
	
	free(nodata);
	
    }
    
    

    finalTime = time(&finalTime);
/*    fclose(fp);
*/

    /* calculate the percentage of valid L2 data points   */
    ll = nscan * npix;
    s = 100.0*(float)pcvalid/(float)ll;
    
    HY1C_out("\nPercent of valid pixels %f\n", s);
    HY1C_out("Accomplished in time = %f [min]\n\n",(float) (finalTime - initialTime)/60.0);
    
 

    if (s < input.vct_fill_pct) {

    	HY1C_out("Failure: %f%% valid L2 ocean data within the L1 scene %s\n\n", s, input.ifile[0]); 

    	if (scan100 <= nscan) {
    	    if ((outtgt = fopen(outtarget,"w")) == NULL) {
             	HY1C_out("%s: Error: Unable to open %s for writing.\n", argv[0], outtarget);
            	exit(NOMATCH_ERROR);
    	    }
    	    fclose(outtgt);
	    sprintf(buf, "rm %s", outtarget);
	    j = strlen(outtarget)+3;
	    buf[j] = '\0';
	    system((const char *)buf); 
	
    	    if ((outar = fopen(outaer,"w")) == NULL) {
           	 HY1C_out("%s: Error: Unable to open %s for writing.\n", argv[0], outaer);
           	 exit(NOMATCH_ERROR);
    	    }
    	    fclose(outar);
	    sprintf(buf, "rm %s", outaer);
	    j = strlen(outaer)+3;
	    buf[j] = '\0';
	    system((const char *)buf); 
	}
		
    	closel1(&l1file);
    	free_l1(&l1rec);
    	free(tgtstr.data);
    	free(aerstr.data);
    	free(tout);
    	free(aout);
	
    	if (TRUE) {
    	    if (TRUE) {	    
	    	free(meta_l3b); 
		free(nbins);
		for (i=0; i<nl2files; i++) {
		     for (j=0; j<nprods; j++) free(sum_prod[i][j]);
		     free(sum_prod[i]);
		     free(sum_chl[i]);
		     free(bins[i]);
		}
		free(sum_prod);
		free(sum_chl);
                free(bins);
		free(numbin);
		free(basebin);
	    } 
	    free(inx);
            free(l2finx);
	    free(l1wavef);
	    free(l2wavef);
       	}
		
	exit(NOMATCH_ERROR);
	
    }
    
    
    
    if (scan100 > nscan) {
    
    	if ((outtgt = fopen(outtarget,"w")) == NULL) {
            HY1C_out("%s: Error: Unable to open %s for writing.\n", argv[0], outtarget);
            exit(FATAL_ERROR);
    	}
    	l2hdr[1] = tglength;
    	fwrite(l2hdr,sizeof(l2hdr),1,outtgt);
    	fseek(outtgt,tglength,SEEK_SET);
    	fwrite((const void *)tout, nscan*tglength, 1, outtgt);
    	fclose(outtgt);
	
    	if ((outar = fopen(outaer,"w")) == NULL) {
            HY1C_out("%s: Error: Unable to open %s for writing.\n", argv[0], outaer);
            exit(FATAL_ERROR);
    	}
    	l2hdr[1] = arlength;
    	fwrite(l2hdr,sizeof(l2hdr),1,outar);
    	fseek(outar,arlength,SEEK_SET);
    	fwrite((const void *)aout, nscan*arlength, 1, outar);
    	fclose(outar);
		
    } else if (oscan % scan100) {

       	if (oscan < OUTSCANLINE) {
 	    if ((outtgt = fopen(outtarget,"w")) == NULL) {
        	HY1C_out("%s: Error: Unable to open %s for writing.\n", argv[0], outtarget);
        	exit(FATAL_ERROR);
    	    }
    	    l2hdr[1] = tglength;
    	    fwrite(l2hdr,sizeof(l2hdr),1,outtgt);
    	    fseek(outtgt,tglength,SEEK_SET);
	} else
	    if ((outtgt = fopen(outtarget,"a")) == NULL) {
        	HY1C_out("%s: Error: Unable to open %s for appending.\n", argv[0], outtarget);
        	exit(FATAL_ERROR);
    	    }
    	fwrite((const void *)tout, wscan*tglength, 1, outtgt);
    	fclose(outtgt);
	
    	if (oscan < OUTSCANLINE) {
    	    if ((outar = fopen(outaer,"w")) == NULL) {
        	HY1C_out("%s: Error: Unable to open %s for writing.\n", argv[0], outaer);
        	exit(FATAL_ERROR);
    	    }
    	    l2hdr[1] = arlength;
    	    fwrite(l2hdr,sizeof(l2hdr),1,outar);
    	    fseek(outar,arlength,SEEK_SET);
	} else
     	    if ((outar = fopen(outaer,"a")) == NULL) {
        	HY1C_out("%s: Error: Unable to open %s for appending.\n", argv[0], outaer);
        	exit(FATAL_ERROR);
	    }
   	fwrite((const void *)aout, wscan*arlength, 1, outar);
    	fclose(outar);
		
    }


    free(tgtstr.data);
    free(aerstr.data);
    free(tout);
    free(aout);
    closel1(&l1file);
    free_l1(&l1rec);
    

    if (TRUE) {
        if (TRUE) {
	    
	    free(meta_l3b); 
	    free(nbins);
	    for (i=0; i<nl2files; i++) {
		for (j=0; j<nprods; j++) free(sum_prod[i][j]);
		free(sum_prod[i]);
		free(sum_chl[i]);
		free(bins[i]);
	    }
	    free(sum_prod);
	    free(sum_chl);
            free(bins);
	    free(numbin);
	    free(basebin);
	} 
	free(inx);
        free(l2finx);
	free(l1wavef);
	free(l2wavef);
    }		


    exit(SUCCESS);
    
    
}








int32 binlocate(int32 *binlist_buf, int32 n_records, int32 bin)
{

   int32  jl = 0, ju = n_records-1, jm;
   int32 i;

   while ( ju - jl > 1 ) {
   	jm = (ju + jl)/2;
	if ((binlist_buf[n_records-1] > binlist_buf[0]) == (bin > binlist_buf[jm]))
	     jl = jm;
	else
	     ju = jm;
   }
   if (jm >= n_records) return( -1 );
   if (binlist_buf[jm] == bin) return(jm);
   if (jm >= 1 && binlist_buf[jm-1] == bin) return(jm-1);
   if (jm+1 < n_records && binlist_buf[jm+1] == bin) return(jm+1);
   return( -1 );
   
}






int32_t diffday(int32_t date1, int32_t date2)
{
  int32_t i;
  int32_t year1, year2;
  int32_t day1, day2, day;

  year1 = date1 / 1000;
  year2 = date2 / 1000;
  day1 = date1 % 1000;
  day2 = date2 % 1000;

  if (year2 > year1) {
     if (isleap(year1) == TRUE) day = 366-day1; else day = 365-day1;
     for (i=year1+1; i<year2; i++)
    	if (isleap(i) == TRUE) day += 366; else day += 365;
     day += day2;
     day = -day;
  } else if (year1 > year2) {
     if (isleap(year2) == TRUE) day = 366-day2; else day = 365-day2;
     for (i=year2+1; i<year1; i++)
    	if (isleap(i) == TRUE) day += 366; else day += 365;
     day += day1;
  } else day = day1 - day2;
  
  return day;

}



void meanday(int32_t year1, int32_t day1, int32_t msec1, int32_t year2, int32_t day2, int32_t msec2, int32_t *myear, int32_t *mday, int32_t *mmsec)
{

  int32_t i;
  int32_t day;
  uint64_t msecs=0, daylong=86400000;


  if (year2 > year1) {
     if (isleap(year1) == TRUE) day = 366-day1; else day = 365-day1;
     for (i=year1+1; i<year2; i++)
    	if (isleap(i) == TRUE) day += 366; else day += 365;
     day += day2;
     msecs = daylong-(uint64_t)msec1;
     for (i=0; i<day-1; i++) msecs += daylong;
     msecs += (uint64_t)msec2;
     msecs /= 2;
     *mmsec = (int32_t)((msec1 + msecs) % daylong);
     *mday = day1 + (int32_t)(msec1 + msecs) / daylong;
     *myear = year1;
     for (i=year1+1; i<=year2; i++) {
	if (isleap(i) == TRUE) if (*mday <= 366) break; else {
	     *mday -= 366; 
	     *myear = i;
	} else if (*mday <= 365) break; else {
	     *mday -= 365; 
	     *myear = i;
	}
     }
  } else if (year1 > year2) {
     if (isleap(year2) == TRUE) day = 366-day2; else day = 365-day2;
     for (i=year2+1; i<year1; i++)
    	if (isleap(i) == TRUE) day += 366; else day += 365;
     day += day1;
     msecs = daylong-(uint64_t)msec2;
     for (i=0; i<day-1; i++) msecs += daylong;
     msecs += (uint64_t)msec1;
     msecs /= 2;
     *mmsec = (int32_t)((msec2 + msecs) % daylong);
     *mday = day2 + (int32_t)(msec2 + msecs) / daylong;
     *myear = year2;
     for (i=year2+1; i<=year1; i++) {
	if (isleap(i) == TRUE) if (*mday <= 366) break; else {
	     *mday -= 366; 
	     *myear = i;
	} else if (*mday <= 365) break; else {
	     *mday -= 365; 
	     *myear = i;
	}
     }
  } else {
     day = labs(day1 - day2);
     if (day1 > day2) {
         msecs = daylong-(uint64_t)msec2;
     	 for (i=0; i<day-1; i++) msecs += daylong;
     	 msecs += (uint64_t)msec1;
	 msecs /= 2;
	 *mmsec = (int32_t)((msec2 + msecs) % daylong);
	 *mday = day2 + (int32_t)(msec2 + msecs) / daylong;
	 *myear = year2;
     }
     else if (day2 > day1) {
     	 msecs = daylong-(uint64_t)msec1;
     	 for (i=0; i<day-1; i++) msecs += daylong;
     	 msecs += (uint64_t)msec2;
	 msecs /= 2;
	 *mmsec = (int32_t)((msec1 + msecs) % daylong);
	 *mday = day1 + (int32_t)(msec1 + msecs) / daylong;
	 *myear = year1;
     } 
     else {
     	 msecs = (uint64_t)labs(msec1 - msec2);
	 msecs /= 2;
	 if (msec1 > msec2) *mmsec = msec2 + msecs; else *mmsec = msec1 + msecs;
	 *mday = day1;
	 *myear = year1;
     }
  }
  
}





uint64_t diffmsec(int32_t year1, int32_t day1, int32_t msec1, int32_t year2, int32_t day2, int32_t msec2)
{

  int32_t i;
  int32_t day;
  uint64_t msecs=0, daylong=86400000;


  if (year2 > year1) {
     if (isleap(year1) == TRUE) day = 366-day1; else day = 365-day1;
     for (i=year1+1; i<year2; i++)
    	if (isleap(i) == TRUE) day += 366; else day += 365;
     day += day2;
     msecs = daylong-(uint64_t)msec1;
     for (i=0; i<day-1; i++) msecs += daylong;
     msecs += (uint64_t)msec2;
  } else if (year1 > year2) {
     if (isleap(year2) == TRUE) day = 366-day2; else day = 365-day2;
     for (i=year2+1; i<year1; i++)
    	if (isleap(i) == TRUE) day += 366; else day += 365;
     day += day1;
     msecs = daylong-(uint64_t)msec2;
     for (i=0; i<day-1; i++) msecs += daylong;
     msecs += (uint64_t)msec1;
  } else {
     day = labs(day1 - day2);
     if (day1 > day2) {
         msecs = daylong-(uint64_t)msec2;
     	 for (i=0; i<day-1; i++) msecs += daylong;
     	 msecs += (uint64_t)msec1;
     }
     else if (day2 > day1) {
     	 msecs = daylong-(uint64_t)msec1;
     	 for (i=0; i<day-1; i++) msecs += daylong;
     	 msecs += (uint64_t)msec2;
     } 
     else msecs = (uint64_t)labs(msec1 - msec2);
  }
     	 
  return( msecs );

}

/* (I'm sorry, but after figuring this routine out, I feel it would be
    good to document what I found for the next traveller, WDR) */

void bandpass_correct(int32_t sensorID1, int32_t sensorID2, float *wave1, float *wave2, int32_t nwave1, int32_t nwave2, float *nLw, float *Lw)
/*******************************************************************

   bandpass_correct

   purpose: interpolate nLw and Lw set from wavelength set of 
            1 instrument to another

   Returns type: void

   Parameters: (in calling order)
      Type              Name            I/O     Description
      ----              ----            ---     -----------
      int32_t           sensorID1        I      incoming sensor ID
      int32_t           sensorID2        I      interpolated sensor ID
      float *           wave1            I      incoming sensor wavelengths
      float *           wave2            I      interpolated sensor wavelengths
      float *           nLw             I/O     incoming radiances, made into
      float *           Lw              I/O     radiances at interpolated 
                                                wavelengths

   Modification history:
      Programmer        Date            Description of change
      ----------        ----            ---------------------
                                        Original development

*******************************************************************/
{
    int32_t iw;

    float nLw_out[NBANDS];
    float  Lw_out[NBANDS];

    float dnLw[NBANDS];
    float dLw [NBANDS];

    
    if (sensorID1 != sensorID2) {
        spline( wave1, nLw, nwave1, 1e30, 1e30, dnLw);
        spline( wave1,  Lw, nwave1, 1e30, 1e30,  dLw);
        for (iw=0; iw<nwave2; iw++) {
  	    splint( wave1, nLw, dnLw, nwave1, wave2[iw], &nLw_out[iw]);
  	    splint( wave1,  Lw, dnLw, nwave1, wave2[iw],  &Lw_out[iw]);
	}
        for (iw=0; iw<nwave2; iw++) {
	    nLw[iw] = MAX(nLw_out[iw],0.0);
	     Lw[iw] = MAX( Lw_out[iw],0.0);
	}
    }        
}


int nLw_SeaWiFS2MODIS(int32_t sensorID, float *nLwbands, float *Lwbands)
{

  static int      firstCall = 1;
  static double   A0[7], A1[7], A2[7], A3[7];
  double  iratio, iratiosq, iratio3, s;
  int    i, nb=NBANDS;
  /* Terra full bandpass from SeaWiFS nominal 11nm bands */
  double  AT0[7]={1.0203581, 0.9774594, 0.9170032, 0.9971813, 0.9056823, 1.0778310, 0.9905358};
  double  AT1[7]={-0.0086241, 0.0077102, 0.0483961, 0.0560854, 0.1494819, -0.1217537, -0.1196560};
  double  AT2[7]={-0.0000673, -0.0033085, -0.0098749, -0.0203777, -0.0156370, 0.0635462, 0.0570984};  
  double  AT3[7]={0.0001394, 0.0004185, 0.0009964, 0.0013208, -0.0001806, -0.0074462, -0.0066622};  
  /* Aqua full bandpass from SeaWiFS nominal 11nm bands  */
  double  AA0[7]={1.0165582, 0.9796995, 0.9213803, 1.0063582, 0.9100613, 1.0755772, 0.9958192};
  double  AA1[7]={-0.0084433, 0.0063936, 0.0436586, 0.0478614, 0.1439513, -0.1206972, -0.1294248};
  double  AA2[7]={-0.0002650, -0.0028901, -0.0091218, -0.0186584, -0.0145033, 0.0631593, 0.0594573};
  double  AA3[7]={0.0001779, 0.0003677, 0.0009270, 0.0011990, -0.0002594, -0.0074087, -0.0068372};
  
  
  /* load parameters */
  if (firstCall) {
  
      if (sensorID == MODIST)
      
      	for (i=0; i<7; i++) { 
	     A0[i] = AT0[i];
      	     A1[i] = AT1[i];
      	     A2[i] = AT2[i];
      	     A3[i] = AT3[i];
        }
	
      else if (sensorID == MODISA)
  
      	for (i=0; i<7; i++) {  
	     A0[i] = AA0[i];
      	     A1[i] = AA1[i];
      	     A2[i] = AA2[i];
      	     A3[i] = AA3[i];
        }
      else {
      
        HY1C_out("nLw_SeaWiFS2MODIS: Error: Needs to be MODISA or MODIST for band pass conversion\n");
        return(1);
      }
	
      firstCall = 0;
       
  }


  if (*(nLwbands+2) <= 0.0 || *(nLwbands+4) <= 0.00001) return(1);
  iratio   = (double) *(nLwbands+2);
  iratio  /= (double) *(nLwbands+4);
  if (iratio > 5.15) return(1);
  iratiosq = iratio*iratio;
  iratio3  = iratiosq*iratio;
  
  for (i=0; i<7; i++) {
  
	s = A0[i]+A1[i]*iratio+A2[i]*iratiosq+A3[i]*iratio3;
	if (i == 3) *(nLwbands+i) = (0.56*(*(nLwbands+i))+0.44*(*(nLwbands+i+1))) * s; else *(nLwbands+i) *= s;
	if (i == 3) *(Lwbands+i)  = (0.56*(*(Lwbands+i)) +0.44*(*(Lwbands+i+1)))  * s; else *(Lwbands+i)  *= s;
	if (*(nLwbands+i) < 0.0) *(nLwbands+i) = 0.0;

  }
  
  if (nb > 9) nb = 9;
  for (i=7; i<nb; i++) {
        if (*(nLwbands+i) < 0.0 && *(nLwbands+i) != -1.0) *(nLwbands+i) = 0.0;
  }
    
  return(0);
    
}  


void tau_interp(int32_t nl1b, float *l1wave, int32_t nl2b, float *l2wave, float *l2tau)
{

  float  *logtau, *tau, angstrom, wgt, *l2wavelog; 
  int    i, ib, i1, i2;
  
  if ( (logtau = (float *)malloc(nl2b*sizeof(float))) == NULL ) {
       HY1C_out("-E- tau_interp: Error allocating memory for L2/L3 file log(AOT)\n");
       exit(FATAL_ERROR);
  }
  if ( (tau = (float *)malloc(nl2b*sizeof(float))) == NULL ) {
       HY1C_out("-E- tau_interp: Error allocating memory for L2/L3 file log(AOT)\n");
       exit(FATAL_ERROR);
  }
  if ( (l2wavelog = (float *)malloc(nl2b*sizeof(float))) == NULL ) {
       HY1C_out("-E- tau_interp: Error allocating memory for L2/L3 file log(wavelength)\n");
       exit(FATAL_ERROR);
  }
  
  for (i=0; i<nl2b; i++) {
      ib = windex(l2wave[i], l1wave, (int)nl1b);
      tau[i] = *(l2tau+ib);
      logtau[i] = (float)log((double)(*(l2tau+ib)));
      l2wavelog[i] = log(l2wave[i]);
  }

  for (i=0; i<nl1b; i++) {
	ib = windex(l1wave[i], l2wave, (int)nl2b);
	if (abs(l1wave[i]-l2wave[ib]) < 2.0) {
	    wgt = 1.0;
	    i2  = ib;
	} else {
	    if (l1wave[i]-l2wave[ib] > 0) {
	       i1 = ib;
	       i2 = ib+1;
	    } else {
	       i1 = ib-1;
	       i2 = ib;
	    }
	    if (i1 < 0) {
	       i1 = 0;
	       i2 = 1;
	    }
	    if (i2 >= nl2b) {
	       i1 = nl2b-2;
	       i2 = nl2b-1;
	    }
	    angstrom = (logtau[i1]-logtau[i2])/(l2wavelog[i1]-l2wavelog[i2]);
	    wgt = pow(l1wave[i]/l2wave[i2], angstrom);
	}
        l2tau[i] = tau[i2] * wgt;
  }
  
  free(logtau);
  free(tau);
  free(l2wavelog);
	      
}




int compar(const void *x, const void *y)
{
   int out;
   uint64_t xx, yy;
   
   memcpy((void *)&xx, (const void *)x, sizeof(uint64_t));
   memcpy((void *)&yy, (const void *)y, sizeof(uint64_t));
   
   if (xx < yy) out=-1; else if (xx > yy) out=1; else out=0;
   
   return( out );
}





int32 open_input_l3b(char *hdf_file, int32 *fid, int32 *sdfid, 
		    int32 *vgid, int32 vdata_id[MAXNVDATA], meta_l3bType *meta_l3b, int32 *nrows)
{
  intn i, j, status;
  int32 vg_ref;
  int32 tag;
  int32 ref;
  int32 vdid;
  char  nam_buf[80];
  char  cls_buf[80];

  if ((*fid = Hopen(hdf_file, DFACC_RDONLY, 0)) < 0)
    {
      HY1C_out("-E- %s line %d: Cannot open input HDF file on Hopen - %s\n", __FILE__,__LINE__, hdf_file);
      return FAIL;
    }

  Vstart(*fid);

  if ((*sdfid = SDstart(hdf_file, DFACC_RDONLY)) < 0)
    {
      HY1C_out("-E- %s line %d: Cannot open input HDF file on SDstart- %s\n", __FILE__,__LINE__, hdf_file);
      return FAIL;
    }

      SDreadattr(*sdfid, SDfindattr(*sdfid, "Start Year"), &meta_l3b->syear);
      SDreadattr(*sdfid, SDfindattr(*sdfid, "Start Day" ), &meta_l3b->sday);
      SDreadattr(*sdfid, SDfindattr(*sdfid, "Start Millisec" ), &meta_l3b->smsec);
      SDreadattr(*sdfid, SDfindattr(*sdfid, "End Year"  ), &meta_l3b->eyear);
      SDreadattr(*sdfid, SDfindattr(*sdfid, "End Day"   ), &meta_l3b->eday);
      SDreadattr(*sdfid, SDfindattr(*sdfid, "End Millisec"   ), &meta_l3b->emsec);
      SDreadattr(*sdfid, SDfindattr(*sdfid, "Sensor Name"), &meta_l3b->sensor_name);

  strcpy(meta_l3b->flag_names, "\x0");
  for (i=0; i<NFLAGS; i++) {
       strcat(meta_l3b->flag_names, l2_flag_lname[i]);
       if (i<NFLAGS-1) strcat(meta_l3b->flag_names, ",\x0");
  }

  *vgid = -1;
  for (i=0; i<MAXNVDATA; i++) {
    vdata_id[i] = -1;
  }
  
  vg_ref = Vfind(*fid, "Level-3 Binned Data"); 
  if (vg_ref > 0) {
    *vgid = Vattach(*fid, vg_ref, "r");

    status = HXsetdir(dirname(hdf_file));
    
    j = 0;
    for (i=0; i<Vntagrefs(*vgid); i++) {
      status = Vgettagref(*vgid, i, &tag, &ref);
      vdid = VSattach(*fid, ref, "r");

      if (vdid == -1) {
	HY1C_out("-E- %s line %d: Problem opening Vdata (reference # %d)\n", __FILE__,__LINE__, ref);
	return FAIL;
      }

      VSgetname(vdid, nam_buf);
      VSgetclass(vdid, cls_buf);


      if (strcmp(cls_buf, "DataMain") == 0) {
	vdata_id[0] = vdid;
      }

      if (strcmp(cls_buf, "Index") == 0) {
	vdata_id[1] = vdid;
	*nrows = VSelts(vdata_id[1]);
      }

      if (strcmp(cls_buf, "DataSubordinate") == 0) {
	vdata_id[2+j++] = vdid;
      }
    }
  }


  return SUCCEED;
}



int32 close_input_l3b(int32 fid, int32 sdfid, int32 vgid, int32 vdata_id[MAXNVDATA])
{

  intn i;

  if (vgid != -1) {
      Vdetach(vgid);

      for (i=0; i<MAXNVDATA; i++) {
	if (vdata_id[i] != -1) VSdetach(vdata_id[i]);
      }
  }

  SDend(sdfid);
  Vend(fid);
  Hclose(fid);


  return SUCCEED;
}

   
   
int empty_tgtaer( int32_t npix, int32_t nscan, int32_t l2sensorID,
  char *outtarget, char *outaer )
/*-----------------------------------------------------------------------------
    empty_trgaer

    Description:  create an empty target and aerosol file

    Arguments:
        Type      Name         I/O   Description
        ----      ----         ---   -----------
        int32_t   npix          I    # pixels to have in file
        int32_t   nscan         I    # lines in the target files
        int32_t   l2sensorID    I    sensor ID of the data
        char *    outtarget     I    target file name
        char *    outaer        I    aerosol file name

    Modification history:

    W. Robinson, SAIC  29 Mar 2011  Original development

----------------------------------------------------------------------------*/
  {
  tgstr tgtstr;
  aestr aerstr;
  int ipx, ibnd, iln;
  int32_t tglength, arlength;
  int32_t l2hdr[5];
  FILE *outtgt, *outar;
 /*
  *  make record lengths for the tgt and aer records
  */
  if( ( tglength = alloc_target( npix, NBANDS, &tgtstr ) ) == 0)
    {
    HY1C_out(
     "%s: Error: Unable to alloc memory for L2 binary target file data %s.\n",
      __FILE__, outtarget);
    exit(FATAL_ERROR);
    }
  if( ( arlength = alloc_aer( npix, NBANDS, &aerstr ) ) == 0)
    {
    HY1C_out(
      "%s: Error: Unable to alloc memory for L2 binary aerosol file data %s.\n",
       __FILE__, outaer);
    exit(FATAL_ERROR);
    }
 /*
  *  set up the header for the files and write them out
  */
  l2hdr[0] = l2sensorID;
  l2hdr[1] = tglength;
  l2hdr[2] = npix;
  l2hdr[3] = 0;
  l2hdr[4] = nscan;

  if( ( outtgt = fopen( outtarget, "w" ) ) == NULL)
    {
    HY1C_out("%s: Error: Unable to open %s for writing.\n", __FILE__, outtarget);
    exit(FATAL_ERROR);
    }
  fwrite( l2hdr, sizeof( l2hdr ), 1, outtgt );
  fseek( outtgt, tglength, SEEK_SET );

  if( ( outar = fopen( outaer, "w" ) ) == NULL )
    {
    HY1C_out("%s: Error: Unable to open %s for writing.\n", __FILE__, outaer);
    exit(FATAL_ERROR);
    }
  l2hdr[1] = arlength;
  fwrite( l2hdr, sizeof(l2hdr), 1, outar );
  fseek( outar, arlength, SEEK_SET );
 /*
  *  zero time info and fill the data values with all -1. or -1
  */
  *( tgtstr.year ) = 0;
  *( tgtstr.day ) = 0;
  *( tgtstr.msec ) = 0;
  for( ipx = 0; ipx < npix; ipx++ )
    {
    for( ibnd = 0; ibnd < NBANDS; ibnd++ )
      {
      *( tgtstr.Lw + ipx + npix + ibnd ) = -1.;
      *( tgtstr.nLw + ipx + npix + ibnd ) = -1.;
      *( aerstr.taua  + ipx + npix + ibnd ) = -1.;
      }
    *( tgtstr.solz + ipx ) = -1.;
    *( aerstr.mod_min + ipx ) = -1;
    *( aerstr.mod_max + ipx ) = -1;
    *( aerstr.mod_rat + ipx ) = -1.;
    }
 /*
  *  write out tgtstr.data aerstr.data to unit
  */
  for( iln = 0; iln < nscan; iln++ )
    {
    fwrite((const void *) tgtstr.data, tglength, 1, outtgt );
    fwrite((const void *) aerstr.data, arlength, 1, outar );
    }
 /*
  *  close the files, free the space, and leave Dodge
  */
  fclose(outtgt);
  fclose(outar);
  free(tgtstr.data);
  free(aerstr.data);

  exit(SUCCESS);
  }
